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Abstract 



Q . We present a Monte Carlo method for the direct evaluation of the differ- 

ence between the free energies of two crystal structures. The method is built 
on a lattice-switch transformation that maps a configuration of one structure 
onto a candidate configuration of the other by 'switching' one set of lattice 
vectors for the other, while keeping the displacements with respect to the lat- 
C^ . tice sites constant. The sampling of the displacement configurations is biased, 

C . multicanonically, to favor paths leading to gateway arrangements for which 

i-^ \ the Monte Carlo switch to the candidate configuration will be accepted. The 

C ' configurations of both structures can then be efficiently sampled in a single 

fi ' process, and the difference between their free energies evaluated from their 

measured probabilities. We explore and exploit the method in the context 
of extensive studies of systems of hard spheres. We show that the efficiency 
JC \ of the method is controlled by the extent to which the switch conserves cor- 

ff^ ' related microstructure. We also show how, microscopically, the procedure 



works: the system finds gateway arrangements which fulfill the sampling bias 
intelligently. We establish, with high precision, the differences between the 
0^ , free energies of the two close packed structures {fee and hep) in both the 

'-"^ \ constant density and the constant pressure ensembles, 

t^ ■ PACS numbers: 05.10.Ln, 65.50.-hm, 64.70.Kb 



I. INTRODUCTION 



Let us pose the problem. We are presented with a material whose chemical composition is 
known; we are provided with a model of the interatomic interactions; and we have identified 
/\ • two candidate crystalline structures. How should we proceed to determine which structure 

j^ ■ will be favored under given conditions? Of course, equilibrium statistical mechanics tells 

us what we must do, in principle: the favored structure will be that which has the greater 
a priori probability or configurational weight; or, in equivalent thermodynamic parlance, 
lower free energy. Thus the task is to compare the configurational weights of (determine the 
difference between the free energies of) the candidate structures. 

A variety of approximate strategies exist for addressing this problem ||l|]. But it is clear 
that if one desires a technique that is both generally applicable and reliable (that is, has 



quantifiable uncertainties) one must look to the Monte Carlo (MC) method , the standard 
computational tool for dealing with many-body systems |^. 

The application of MC methods to the study of phase-behavior presents a generic prob- 
lem @J^: the free energy of a phase cannot be expressed (in a practically useful form) as a 
canonical average over the associated configurations; free-energy-estimation inevitably en- 
tails simulations that visit a substantially wider spectrum of configurations, which together 
form a path through configuration space 0. The strategic choices to be made concern the 
path itself — ultimately, the physical character of the additional configurations sampled — 
and the sampling procedure. 

An acceptable path will fall into one or other of two categories — we will call them 
reference- state and inter-phase paths. A reference-state path links (comprises sets of config- 
urations that interpolate between) the configuration space associated with each phase to the 
configuration space associated with some reference system whose free energy is known. 
An inter-phase path links the configuration space of one phase to that of the other. Both 
categories of path embrace many further sub-categories. Thus, a reference-state-path may 
run through a space of thermodynamic coordinates or through a space of model parameters. 
An inter-phase path may be 'physically-motivated' — modeling authentically the configura- 
tions through which a system actually passes in the course of a phase transformation; or it 
may be 'computationally-motivated' ('non-physical') — designed, pragmatically, to deliver a 
result. 

The sampling procedures used to explore the chosen path also fall, broadly, into one or 
other of two categories - we will call them multi-stage and single-stage. The multi-stage 
approach entails a series of independent simulations each of which explores a different point 
on the path; the simulations may determine simply the derivative of the free energy at each 
point (the integration method, IM), or the difference between the free energies of adjacent 
points (the overlap method). The single-stage approach involves, in essence one simulation 
exploring the entire path. 

There are very many ways in which one can respond to these strategic choices. Many of 



them are represented in the large literature devoted to this problem [§- [|14[. But all of them, 
in our view, lack one or more of the characteristics (generality, transparency, precision) of a 
definitively-satisfactory solution to such a fundamental and simply-posed problem. 

In seeking that solution it seems to us there are good a priori grounds for favoring an 
inter-phase path, explored by single-stage sampling. The prejudice on the choice of path 
reflects the fact that, in using a reference state path, one has to determine, separately, 
the absolute free energies of each phase. These absolute free energies are typically very 
large -arbitrarily so in the vicinity of a phase boundary- compared to the quantity (their 
difference) which is actually of interest. In contrast, using an inter-phase path allows one 
to focus directly on this quantity. The a priori preference for a single-stage sampling rests 
on the transparency with which the associated uncertainties (error bounds) are prescribed. 
We shall return to these points in section 0. With these strategic choices made, one is left 
with two tasks — one conceptual (designing the inter-phase path), and the other practical 
(formulating the sampling algorithm). 

The practical issue is relatively easily addressed. In recent years, the Monte Carlo toolkit 
has been significantly enhanced to provide a range of extended sampling techniques — multi- 



canonical [0, expanded-ensemble []TB[ and simulated-tempering |T^. These methods (whose 



origins can be traced back to much earlier pioneering work [0) allow one to construct a 
MC procedure that will traverse virtually any desired path through configuration space. 
Here we adopt the multicanonical framework. In this framework, the desired path is rep- 
resented as a discrete series of macrostates, defined by some chosen macroscopic property 
0]; in the multicanonical sampling procedure each macrostate is visited with a probability 
that is enhanced, or diminished, with respect to its canonical value, by an amount that 
is controlled by a multicanonical weight; the set of weights is constructed so that, while 
the canonical probabilities vary vastly over the path, the multicanonical probabilities are 
essentially constant, allowing the whole path to be negotiated in one simulation. 

The core issue is, then, the design of the inter-phase-path — at heart, the choice of an 
appropriate order parameter [|l^. The choice is important: it determines, implicitly |2^, the 
nature of the configurations sampled during the inter-phase traverse, the MC-time required 
for that traverse, and thence the statistical quality of the final results. 

Outside the context of structural-phase behavior -in the case of liquid-gas phase behavior 
for example — the choice is clear and a multicanonical strategy is securely in place. The 
order-parameter is identified with that — the density — associated with the accompanying 
critical point. The resulting inter-phase configurations are then generically inhomogeneous, 
comprising two coexisting regions (one of each phase), separated by an interface. On this 
path, it is the free energy cost of this interface that provides the ergodic barrier which has 
to be surmounted by multicanonical weighting ||2l]] . The passage along the path (the motion 
of the interface) involves processes which differ only in scale from those already represented 
in the microscopic dynamics of a single phase. This approach is illustrated schematically in 
Fig. (|I]a). It has been successfully used in studies of phase behavior in ferromagnets p2| . 



fluids [O and lattice gauge theories [24 



In the context of structural phase behavior it is clear that this kind of strategy will not 
usually be fruitful [^. In such systems a traverse through an inhomogeneous two-phase 
(necessarily non-crystalline) region will involve substantial, physically slow, restructuring 
— vulnerable to further ergodic traps, and compounding the intrinsic slowness of the mul- 
ticanonical sampling process. To the two general a priori preferences expressed above we 
thus add a third, speciflc to the structural context: the inter-phase path should comprise 
macrostates that are single-phase, and crystalline. This paper shows how to identify, build 
and exploit a path of this kind. 

The key ideas are simple. In any crystalline conflguration each atomic position coordinate 
may be expressed as the sum of a lattice vector and a displacement vector. The conflgura- 
tion space associated with each structure, individually, may be explored by standard MC 
procedures which stochastically update the displacement vectors while keeping the lattice 
vectors constant. In principle the passage from one phase to the other may be accomplished 
by a lattice switch (LS) in which one entire set of lattice vectors is replaced with the other, 
while the displacement vectors are held flxed. Formally this LS can be incorporated into the 
MC procedure simply by treating the lattice type as an additional stochastic variable. In 
practice this inter-phase 'path' (blind-leap) will not work. Implemented this way the LS will 
map a 'typical' conflguration of one structure onto an 'untypical' (high-energy) conjugate 
conflguration [^; the associated MC step will generally be rejected. To make it work the 
LS needs to be extended to include two segments of 'path' (each lying entirely within one 
phase) which connect the sets of equilibrium conflgurations with the special conflgurations 



(we will call them gateway configurations) from which a successful LS can be initiated p8 . 
These path-segments may be labeled by an 'order-parameter' which measures the mismatch 
between the energies of the configurations linked by LS. This order parameter has a high 
value for the equilibrium configurations, lying at one end of a path segment: these configu- 
rations are not energy-matched [^ to their conjugates. It has a low value for the gateway 
configurations at the other end: gateway configurations (whatever other attributes they 
may have) are necessarily energy-matched to their conjugates. Multicanonical weights are 
attached to the macrostates of this order parameter, so that the multicanonical sampling 
procedure explores both path segments evenly, surmounting the probabilistic barrier which, 
in this case, reflects the smallness of the statistical weight of the gateway configurations. 
Together, the multicanonical sampling and the lattice switch provide a configuration-space 
'look and leap' (Fig. |I]b) which visits both phases while remaining at all times crystalline. 

The LS method was introduced by us and described in outline form in an earlier brief 
communication |]3D|. Since that time it has been applied by two other groups [|T^^. The 
present paper has three principal objectives. 

The first objective (with which we have already engaged in the preceding discussion) is to 
explain the core idea more fully: the 'idea' (biased sampling to facilitate a global coordinate 
change) represents, we believe, a significantly new form of extended sampling, which merits 
further exposure. 

Our second objective is to achieve a deeper understanding of how the process works -in 
particular the implications of the form chosen for the LS operation adopted (it is not unique) 
and the microscopic character of the gateway configurations which the system locates in 
response to the multicanonical weighting, tailored to support that operation. We show that 
the efficiency of the LS operation depends significantly on the extent to which it conserves 
correlated microstructure. And we find that the gateway configurations have features which 
reflect the specific nature of the lattice-switch transformation we adopt, in a microscopically 
intelligible (even intelligent) way. 

Our third objective is to extend our study of the phase-behavior of hard-spheres. This 
problem is of enduring interest, displaying a richness that belies the simplicity of the model 
itself ||3^. The relative stability of the two closed-packed {fee and hep) structures is par- 
ticularly finely balanced: the entropy difference [^] is so small (smaller than the entropy 



change at freezing by of order 10 ^) that it can easily be lost in statistical uncertainties. 
Discrepancies (large, in relative terms) between a recent IM study [|10| and its predeces- 



sors provided the motivation for our development of the LS method. In this paper we 
present results in the constant-density ensemble, both near the melting density and at the 
close-packed limit. In so doing we resolve the discrepancy between the results, near melting, 
reported in our initial study [Q and those — also using LS — reported recently by Pronk & 
Frenkel |^l|: the fault was ours, stemming from a failure to recognize the consequences of 



center-of-mass drift. We also show that the method can be extended straightforwardly -in 
this case at least- to the constant-pressure ensemble. 

The paper is structured as follows. Section y sets out the theoretical framework. We 
define the model, the competing structures, and the associated configurational weights: in 
the case of hard sphere systems the latter are purely entropic. We identify an appropriate 
form of lattice-switch transformation: here, it is designed to capitalize on the close-packed 
layers common to both structures. To bias the displacement sampling we need to define an 



appropriate measure of the 'energy cost' of the lattice switch; we will see that the number 
of pairs of overlapping spheres created by the transformation fulfills this role simply and 
effectively. The efficiency of the method also potentially depends on the choice of represen- 
tation of both the lattice-to-lattice mapping and the particle displacements: we discuss the 
principles involved in the choice of representation. Section ^ provides computational im- 
plementation details, including the procedures used to evolve an appropriate multicanonical 
sampling distribution. Section ^V\ contains our results. Finally, in section |V|, we offer our 
conclusions in relation to both the hard sphere system and the lattice-switch method. 

II. FORMULATION 

A. The model system 

We consider a system of N particles, of spatial coordinates {r}, confined within a volume 
V, and subject to periodic boundary conditions. The interactions are those of hard spheres 
of diameter D] the configurational energy is of the form 

^^ ^^ 1 c>o otherwise ^ ^ 

where Vij =\ fi — fj \. The total configurational weight associated with this system is 

n{N,V) = l[[[ dfil n Q{r^,-D) (2) 

where G(x) = 1 (0) for x > (< 0), and the product on < ij > extends over all particle 
pairs. The associated entropy density is 

s{N,V)^^\nn{N,V) (3) 

We are concerned with the entropy of specific phases (the two familiar crystalline close- 
packed structures) of this system. In general, the entropy of a phase measures the weight 
of the configurations satisfying some constraint that is characteristic of that phase. It 
is therefore necessary in principle (although in practice the issue is typically skirted) to 
formulate a constraint that identifies a configuration as 'belonging to' a given crystalline 
phase. One can do so -very naturally, and in the traditions of lattice dynamics ||3^- by 
decomposing the particle position coordinates into a sum of 'lattice' and 'displacement' 
vectors: 

fi = R^ + Ui (4) 

Here {R}a = -R", i = 1 ... iV is a set of fixed (configuration-independent) vectors associated 
with the crystalline structure labeled a. We will refer to them as 'lattice vectors'. But 
we use this term a little loosely: more precisely, we mean the set of vectors identified by 
the orthodox crystallographic lattice, convolved with the orthodox basis |^. The other 
vectors, {u} , represent displacements with respect to the 'lattice' sites; the symmetry of 



the structures of interest here ensures that these displacements have zero ensemble average. 
This framework provides us with a number of ways of identifying the configurations to be 
associated with structure a. First, one might adopt the criterion that all particle displace- 
ments, with respect to the associated lattice sites, lie within some nominated spatial cutoff, 
chosen to be sufficiently large that the results are independent of its specific value. This 
criterion has the merit that it does not stray beyond the concepts of equilibrium statistical 
mechanics. Alternatively one might identify the relevant configurations as the set that are 
accessible from some member of the set (the perfect crystalline state, for example) within 
some nominated temporal cutoff. The merit of this choice is that it is a quasi-formal expres- 
sion of what, in practice, computer simulation attacks on this problem actually do, albeit 
implicitly: the free energy assigned to a phase (in, for example, IM-based studies) repre- 
sents the weight of the configuration space sampled on the time scale of the simulation. The 
result should be independent of that time scale provided it (the scale) is long enough that 
the configuration space of each structure is effectively sampled, but still short compared to 
inter-phase crossing times. Whichever view one takes (in practice we adopt the latter: see 
section p.11 A] ) one may write, for the configurational weight associated with structure a 



niN, V, a) = n[ L^^^] n 0(^^. - D] 



(5) 



<«J> 



where /rrn signifies integration subject to the chosen configurational constraint. 



In the thermodynamic (A^ —> oo) limit, the associated entropy density 

siN,V,a)^^Ann{N,V,a) 



N 
depends only on the particle number density, which we write in the dimensionless form 



(6) 



(7) 



Pep V2/D3 
where pcp, the number density at close packing, provides the natural scale. The range of 



interest to us here extends from the melting density p ~ 0.736 ||3^ through to the close- 
packed limit p = 1. 

In the close-packed limit the configurational integral (Eq. ^ may be rewritten [^ as the 
product of two terms: 



n{N,v,a) = noiN,v)nc 

The first term here is defined by 



no{N,v) 



■ De - 
l-e. 



3N 



(8) 



(9a) 



with 



= 1 - pi/3 



(9b) 



The associated contribution to the entropy is logarithmically divergent in the close- 



packed limit ^9[, but independent of the phase. The second contribution to the configura- 



tional integral is defined by 

„ nn 

Qa = IlijL^^^] n ^(4 + 1) [1 + 0(e)] (10a) 



<»j> 



where |^| 



Uij =Ui — Uj = Uijnfj + uj-j (10b) 

while n°'j is a unit vector from lattice site j to nearest neighbor lattice site i. The associated 
contribution to the entropy is finite, but depends on the phase through the geometry of the 



nearest neighbor vectors. It may be visualized as that of a set of hard dodecahedra ||41|| . 

Now let us recall that the quantity of immediate interest is the difference between the 
entropy-densities of the two phases. It may be written as 

Aso^p = s{N,V,a) - s{N,V,P) = ^\n'JZ^p{N,V) (11) 

where 

^"^^^' ^^ = niN,V,(3) - P(/3 I AT, V) ^^^^ 

Here P{a \ N, V) is the probability that a system, free to explore the joint configuration 
space of the two structures (and visiting configurations with the appropriate probabilities 
— all equal in this case) will be found in some configuration of structure a. 

In the constant density ensemble, then, the computational task is to determine the ratio 
defined by Eq. (|T^. In the constant pressure ensemble we require the ratio 7laf3{N,P*) of 
the partition functions 

Z{N, P\ a)= I dVn{N, V, a)e-^*^ (13) 



where P* is a measure of the pressure |^^. The associated thermodynamic potential is the 
gibbs free energy density defined by 

g{N, P\ a) = -^ In Z{N, P*, a) (14) 



so that, in analogy with Eq. ([TT|), 



Ag^^ ^ g{N, P*, a) - g{N, P*, P) = -1 In 7^„^(iV, P*) (15) 



B. The lattice-switch method 

The two close-packed structures of interest here are shown schematically in Fig. (^. 
In principle there are many transformations which will map one set of lattice vectors into 
the other; we shall consider the criteria guiding the choice in section |IIC|. The mapping 



used in most of the work reported here is shown schematically in Fig. (^). This scheme 
exploits the fact that the two structures differ only in respect of the stacking pattern of 
the close-packed planes. A suitable transformation can then be constructed that entails, 
simply, translating appropriate close-packed planes. By 'translate' we mean, more precisely, 
'relocate at a position defined by an appropriate translation vector': one should not think 
of the planes as 'sliding through' the intermediate states. 

Figure (|]) shows the application only to the ]5er/ect-crystal-configurations where energy- 
matching is guaranteed |^. In general (that is, for 'typical' configurations: see Fig. (^ 
for an example) the two configurations related by the LS operation will not be energy- 
matched: since adjacent planes are translated differently, the translations may — indeed, 
with overwhelming probability will — map a realizable configuration (of one structure), in 
which there are no overlapping spheres, onto an unrealizable configuration (of the other) in 
which there are overlaps. A MC lattice switch 'move' will be rejected for most configurations. 



But not quite all: gateway configurations (configurations that are energy-matched p9| to 
their conjugates) must exist, in significant measure. In particular, it is clear on grounds 
of continuity that configurations 'close enough' to perfect-crystal form must fall into this 
category. One might therefore choose these 'small-displacement' configurations to act as 
the gateway states, and define a multicanonical weighting procedure accordingly. However, 
one can avoid having to make this explicit choice, and, instead, let the system find gateway 
configurations itself. To do so we must define a measure of the mismatch between the 
energies of the configurations linked by the transformation. 

In the present context that mismatch is quantified by the number of pairs of overlapping 
spheres created by the transformation. To that end let M{{u},a) denote the number of 
overlapping pairs associated with the displacements {u} within the structure a. And define 

a 



Mi{u}) = M({n}, hep) - M{{u}Jee) (16) 

Since M{{u}, a) will necessarily be zero for any realizable set of displacements of structure 
a, the overlap order parameter Ai is necessarily > (< 0) for realizable configurations of 
the fee (hep) structure. Figure (^ provides a concrete example. The gateway configura- 
tions may then be identified {without prejudging their microscopic character) as the set of 
configurations for which A^ = 0: a displacement pattern {u} for which A^ = is realizable 
in both structures (energy- matched). A LS MC step initiated from an A^ = configuration 
will be accepted; if initiated from outside this set of configurations it will be rejected. 

The sampling algorithm must thus be multicanonically customized so as to enhance the 
probability along a notional line in A^-space, extending from the 'equilibrium' Al-values 
(reflecting the number of overlaps created by a LS acting on a typieal configuration) through 
to the Al = gateway configurations. This aim is realized by augmenting the system energy 
function Eq. (|T]): 

E{{f}) -. E{{r}) + r,{M{{u})) ^ E{{r}) (17) 



where ri{M.),M. = 0, ±1, ±2 . . . constitute a set of multicanonical weights |]T5|. These 
weights need to be chosen so as to allow the system to access the 7W = gateway con- 
figurations, and thence (through the LS) the full joint configuration space of the two struc- 
tures. The desired ratio of configurational weights, which refiects the canonical distribution 
P{Ai I A^, V) (Eq. |T2|) may then be estimated from the measured multicanonical distiihution, 
P{M I A^, V, {r]{M)]) with the identification 

V (-r in - ^M>oPiM I N,V) _ EM>oPiM I N,V,{viM)})e^^^^ 

^^'^'"^"^ ' ^ j:M<oPiM I N,V) EM<oPiM I N,V,{viM)})e^(M) ^'^) 

Here the exponential re-weighting of the multicanonical distribution folds out the bias asso- 
ciated with the weights, whose residual effects are then simply as desired — the removal of 
the ergodic barrier between the two branches of the distribution. 



C. Representations: tuning the lattice switch 

We have presented the LS method in its simplest realization -the one we have used for 
most of the studies reported here. We now outline two important respects [^ in which 



some degree of generalization is possible, and may be desirable, in subsequent applications. 
Both involve the choice of representation of the LS transformation. 

We have already alluded to the first point: there are many forms of lattice-to-lattice 
mapping. It is clear that the efficiency of the method will depend significantly upon the 
mapping chosen. Evidently the choice should be made so as to match up, as closely as 
possible, the energy of the two configurations it links. In the context of hard spheres this 
aim is realized by choosing the mapping which gives the smallest equilibrium overlap-count 
(mean | Ai |-value), which gives a measure of the entropic barrier that has to be negotiated 
by the multicanonical procedure. The smaller this barrier, the shorter is the path to the 
gateway configurations from which a successful LS may be launched. Since the multicanon- 
ical simulations traverse this path only slowly (essentially diffusively, at best) the gains here 
are potentially substantial. It is intuitively clear that the scheme described above will fulfill 
this criterion well: in this representation, the LS translates close-packed planes bodily, so it 
can create overlaps only between spheres associated with different planes. But it is useful to 
explore other schemes -partly to check that there is no significantly better alternative, but 
principally to understand the different factors that control the efficiency. We have done so; 



the results are to be found in section IV A 



There is a second - less obvious- generalization of the framework. In the simple real- 
ization, the particle positions are written in the 'lattice plus displacement' representation 
provided by Eq. (^. The LS operation then maps a configuration of one structure onto a 
configuration of the other with the same set of displacements. This is unnecessarily restric- 
tive. More generally we are at liberty to write, in place of Eq. (^), 

f=^" + T"-M (19) 

where r, i?° and u are now column vectors with 3A^ elements and T° is a 3A^ x 3A^ non- 
singular matrix, whose form (possibly {u}-dependent) is at our disposal. Eq. (^) is then 
replaced by 

9 



n{N, V, a) = Uljr^du,] ■ det T" H ^i^ij - D) (20) 



"' <^3> 



From the standpoint of the (standard) single-phase part of the MC procedure, this change 
in representation is equivalent to changing the form of the configurational energy: 

E{{r}) ^ E{{r}) - In [det T'^] (21) 

This change introduces some computational overheads, which could be substantial if the T- 
transformation is not local. The potential pay-off lies in the LS part of the MC procedure. 
One might hope to be able to tune the form of the T-matrix so that 'typical' configurations 
of the one structure are mapped (by LS) into 'typical' configurations of the other. In the 
case of the hard sphere problem, however, our results (section |1VA]) suggest that there is 
little to be gained here by this kind of tuning. 



III. IMPLEMENTATION 

A. Monte Carlo procedures 

First we consider the procedure for MC-sampling of the particle displacements, for a 
given structure (set of lattice vectors). As discussed in section [II 7\] this sampling should. 



in principle, satisfy some appropriate configurational constraint [Q. In our original studies 



|30| we chose to implement this constraint explicitly^ through our sampling distribution: 



candidate displacements were drawn from a flat ('top-hat') distribution. This procedure can 
be made to work. But the constraint explicitly breaks the translational invariance; and one 
must deal with the consequences. In particular the configurational integral effectively being 
evaluated then depends upon the location of the center of mass and thence upon the top-hat 
cut-off; this dependence sets in when the displacement acquired by the center of mass, in 
the course of its slow diffusive motion, becomes comparable with the top-hat-cut-off. One 
can avoid this problem simply by fixing the center of mass. Our failure to do so in |3^ led 



to results which differ significantly from those we present here. In the studies reported here 
we have chosen the 'implicit' realization of the configurational constraint (practically, but 
not conceptually equivalent to ignoring it) which rests (section [11 A|) on time scales. Spheres 
were chosen at random, and ifisX- changes io the current displacement drawn from a uniform 
distribution. The displacement update is accepted according to the Metropolis prescription 



Pa{{u} -^ {u'}) = min {l, exp [-AE{{f})] } (22) 

where E{{f}) is defined in Eq. (|1^. In addition to the constraint that the update should 
yield a realizable configuration of the current phase, this acceptance probability reflects the 
chosen weights which are defined (section [1IIB| explains how) on the space of the overlap 
order parameter Ai (Eq. |l^). To minimize the computational time spent determining how 
a proposed move affects the value of Ai we used a local overlap array, holding information 
on which neighbors of a given sphere currently overlap with that sphere in the conjugate 
configuration generated by a LS. 

10 



The representation of the close-packed hmit provided by Eq. (|10a]) can be handled with 
only minor amendments: the constraint Vij > D identifying realizable configurations is re- 
placed by a constraint on the scaled displacement-difference coordinates, m|'- > — 1. The over- 
lap order parameter (measuring the number of times the hard sphere constraint is violated 
in the conjugate configuration) is redefined accordingly. In this limit particle 'interactions' 
(encounters) may occur only between immediate neighbors. At other densities we allowed 
for the possibility of encounters between nominal second-neighbors. We found however that 
although the number of such encounters grows rapidly with the approach to the melting 
density, the consequences for the relative entropy of the two structures is insignificant under 
the conditions studied here [^ . 



In addition to particle moves the constant-pressure simulations require updates of the 
simulation-cell-parameters. In such an update (implemented on average once per sweep) a 
trial set of cell-parameters are selected, and accepted with probability 



Pa{V -^ V) = min {l, exp [-AE{{f}) - P*AV + N \n{V' /V)] } (23) 

where V is the volume associated with the trial parameters. Note that this kind of update 
-a dilation- changes E{{f}) both trivially (so as to forbid moves causing 'real' overlaps) and 
more subtly through changes in the count of the overlaps in the conjugate configuration. A 
volume update thus requires recalculation of the entire local overlap array. 

Now consider the lattice switch. The switch may be viewed as an updating of the 
'lattice '-type a, regarded as a stochastic variable. The prescription for such an update is 
quite simple. After every particle update the value of Ai is checked (it is already known). 
The LS is performed if (and only if) the gateway condition A^ = is satisfied. 

B. Calculating the weights 



The determination of an acceptable set of multicanonical weights |15| can be accom- 
plished in a number of ways — none, seemingly, entirely systematic. We describe briefly the 
techniques we have used in the present study. Figure (^) provides some illustration. For 
further details and references to other work the reader is referred to [p!5| , p6| , |50| , pT| . 

The simplest method is the visited states (VS) technique 0]. In this approach a suitable 
set of weights is evolved through an iterative process (Fig. ^, the next set of weights de- 
pending upon the distribution of the (overlap) order parameter over the macrostates visited 
using the current set of weights. This process is repeated until the weights yield a distribu- 
tion P{A4 I N,V,{7]{Ai)}) that is effectively flat. This method proved quite adequate for 
our smallest system. 

For larger systems, however, we found it more efficient to appeal to the transition proba- 
bility (TP) method |^. In the simplest realization of this method the simulation is initiated 



from a 'cold' (zero-displacement) configuration (a member of the A^ = macrostate) for one 
structure. In the course of its subsequent evolution towards equilibrium for that structure 
the numbers of transitions between different A^-macrostates are recorded, and subsequently 
used to construct an estimator of the macrostate-transition-probability matrix. This TP 
matrix can be used to estimate the macrostate probability distribution and thence to pro- 
vide an estimate for a set of weights (Fig. ^, which can in turn be refined via VS. For our 
intermediate size system this method worked well. 
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In the case of our largest system we found it necessary to modify the method somewhat, 
so as to hmit the rate at which the simulation passes through A^-space. One way of doing 
this is to constrain the system to macrostates with overlap order parameters below some 
'barrier'- value, which is gradually incremented (moved 'out' in A^-space), at intervals of the 
order of the equilibration time. 

By fiat the two structures have the same weights for A^ = 0. In principle, the weights 
associated with the two structures for non-zero | Al | are different (i.e. r]{J\4) ^ r]{—M)^, 
and have to be evolved separately. In practice the weights of the two structures are very 
similar — a reflection of the similarity of the entropies of the two phases. Consequently, one 
set of weights provides an excellent first approximation to the other, for refinement by VS. 

C. Simulation details 

The specific form of the LS operation we have chosen (Fig. 0) imposes restrictions on 
the geometry of the system simulated: with normal periodic boundary conditions the sys- 
tem must comprise integral multiples of 6 close-packed planes. It is possible to avoid this 
restriction by using more elaborate boundary conditions [1^ , but we chose to avoid this com- 



plication and simulate systems comprising 6^ = 216, 12^ = 1728 and 18^ = 5832 spheres. 
Simulations were performed at two densities, namely (cf Eq. 0) p = 0.7778 [^, and p = 1, 
the close-packed limit. 

The maximum step size for displacement updating was chosen so as to minimize the 
autocorrelation time of the overlap order parameter (Eq. |16D. We found a maximum step 
size of O.ISD produced the best results at p = 0.7778, while a value close to unity was found 
to be appropriate in the close-packed limit, in the representation (and scaled units) given in 
Eq. ([lOil). 



A significant proportion of our simulation time was devoted to the process of weight- 
determination. For our largest system we used 10^ Monte Carlo sweeps (MCS) to generate a 
first (TP) estimate of the weights, with a further 5 x 10^ MCS devoted to weight-refinement 
using VS. 

The free-energy differences of interest were then determined by further simulations in 
the multicanonically-weighted ensemble. For each system (density, and size) we performed 
a series of runs each long on the scale of the autocorrelation time of the overlap order 
parameter. Each of these runs then provides an independent estimate of the (logarithm of 
the) probability ratio required (Eq.|TB|). The standard deviation of these estimates provides 
a basis for assigning an associated statistical uncertainty. Implementing this stage required 
simulation times ranging from ~ 2.5 x 10^ MCS for A^ = 216 to ~ 4 x 10^ MCS for A^ = 5832. 

IV. RESULTS 
A. The effects of the representation 

As discussed in section |IIC| the LS operation can be implemented with different choices 



of representation of the lattice-mapping or the particle displacements [Eq 
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The efficiency of a lattice mapping is measured (inversely) by the equilibrium overlap 
count. Table H shows results for a variety of mappings, chosen to expose the different 
factors that control the mapping efficiency. Mapping number 1 is the one described in 
Fig. I, and used throughout this work: the notation (0, — t, +i) signifies that the three pairs 
of planes counting from the top of Fig. (^ are translated respectively by 0, —t and +t. 
A similar convention is used to label mappings 2 and 3. In mapping 4 ('random-plane') 
an hep configuration is generated by taking an fee configuration and restacking its close- 
packed planes in a random order, in an hep pattern. In mapping 5 ('random-site') an hep 
configuration is generated by mapping the particle displacements in an fee configuration 
randomly on to the sites of an hep lattice. 

The random-site mapping (number 5) shows the largest overlap count. One can account 
for its value, rather well, by regarding the particle displacements as isotropie, gaussian and 
independent of structure [^], and estimating the probability that two particles associated 
with nearest-neighbor sites, and with displacements drawn randomly from this distribution, 
will overlap. 

Using the random-plane mapping (number 4) cuts the overlap count by a factor of (a 
little more than) two with respect to random-site. This efficiency gain simply reflects the 
fact that of the 6A^ potential overlaps between near-neighbors, only the 3A^ associated with 
neighbors in different (but adjacent) planes can now contribute. 

Mapping 3 simply generates one fee configuration from another (it is useful only because 
it is informative): its overlap count is cut by a further factor of two. This reflects the 
fact that this mapping (like mappings 1 and 2) moves close-packed planes in pairs, thus 
guaranteeing no overlaps between the two members of each pair. 

Mappings 2 and 1 show further — smaller but still practically useful — cuts in the overlap 
count. The origin of these gains is more interesting. It is clear that they must reflect the 
size of the translation vector used: mappings 1 through 3 differ only in this respect. This 
vector controls the extent of the shear which the mapping introduces between successive 
pairs of planes. The following interpretation seems reasonable. The displacement patterns 
in adjacent planes will be correlated to some extent, with undulations in one surface (the 
z-components of the displacements) matched to undulations in its neighbor. The smaller 
the shear, the more closely these undulations will remain matched to one another (in the 
conjugate configuration), and the smaller the overlap count. With increasing shear, this 
advantage is lost and the behavior should (and indeed does) approach the limit (one quarter 
of the overlap count for mapping 5) one would expect in the absence of such correlations. The 
fact that this 'approach' is already apparent in the performance of mapping 2 is consistent 
with the fact that the measured correlation length of the surface undulations at the density 
concerned is found to be close to the magnitude of the translation vector t. 

These results help to clarify the factors which control the overlap count of the mapping 
(number 1) we have actually used. It is tempting to attribute the overlaps to the fact that 
the LS {fee — * hep, say) maps each particle from an environment in which adjacent close- 
packed planes have different stacking labels (A and C, say) to one in which they have the 
same label (C, say). The results for mappings 1-3 show that it would be misleading to think 
this way. The overlaps simply reflect the numbers of particles that 'see' a new adjacent 
close packed plane (irrespective of its label), and the extent to which it is 'new'. This is the 
reason for the similarity between the overlap counts for the two structures (section |1V C| ) . It 
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shows, moreover, that any simple [^ tuning of the displacement representation (the choice 



of T- matrix) is likely to be of no advantage here [Ml - 

B. How it works: the gateway configurations 

A LS operation will work (be accepted as a MC move) only when launched from a small 
subset of the configurations actually visited: these, by definition, are the 'gateway config- 
urations'. As noted earlier, one could identify a priori configurations (those characterized 
by 'small enough' displacements) which fall into this set. But we have elected, rather, to 
let the system (the algorithm) identify them on the basis of their defining characteristic — 
that they have zero overlap order parameter M. ||2^ . It is then interesting to investigate the 
microscopic characteristics of the configurations picked out by this constraint. Figure (^ 
shows the distribution of the separation, d, between adjacent close-packed (x — y) planes |]55| , 
for A^-macrostates corresponding to equilibrium fee, equihbrium hep, and gateway {A4 = 0) 
regions. The macrostates corresponding to the equilibrium crystal structures have similar, 
near-gaussian, ^-distributions. In contrast, for the gateway macrostate the distribution is 
bi-modal: in this macrostate, some planes are systematically moved closer to one another, 
while (in equal measure) others are shifted apart. On closer examination one finds that it 
is the planes which are translated together by the LS (eg the pair of planes marked (i) in 
Fig. (^) that fall into the first category, while the planes that are translated differently hy 
the LS (eg the pair of planes marked (ii) in Fig. (^)) fall into the second. The evolution, 
with M. , of the mean plane separation (for both categories) is shown in Fig. d^). The 
behavior thus unearthed is entirely reasonable. The LS operation can only create overlaps 
between neighboring planes which are translated by different amounts (sheared with respect 
to one another). The algorithm resolves the task set by the bias towards 7W = by moving 
these pairs of planes (the ones vulnerable to overlaps) further apart, at the expense of a 
compression of the others ||56|. In simulations conducted at constant pressure this effect 



(still present) is supported by a second. Fig. (^) shows that the algorithm now exploits 
the additional degrees of freedom (the shape of the simulation cell) to locate gateway states 
with values of the c/a ratio enhanced above the ideal close-packed value |^8[. Again, the 
advantages with respect to the switch are clear. 

It is tempting to say that the sampling is intelligent. In any event it is clear that the 
algorithm locates and utilizes configurations which it would be difficult to exploit explicitly 
in the design of the switch operation. 

C. Entropies of crystalline structures 

The essential output of a LS-simulation is in the form of the normalized probability distri- 
bution of the overlap-order-parameter, reweighted to remove the bias in the multicanonically- 
weighted distribution actually measured. Figure (H) shows the results for this distribution 
(at p = 0.7778) for three different N values. As one would expect the distributions each 
comprise two peaks (one associated with each phase) each of which is nearly gaussian |5^ 



and sharpens with increasing A^ |B^. Note the close correspondence between the equilib- 
rium overlap counts for the two structures. This result is not required by definition, or any 
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obvious symmetry. Rather it should be seen as a further manifestation (the smallness of the 
entropy difference between the phases is the prime one) of the similarity of the local particle 
environments in the two structures. 

The relative weights of the two peaks is a direct measure of the difference between the 



entropies of the two structures (Eqs. ITT], O, 11^). Since the entropies are extensive the ratio 



of the peak weights grows exponentially with A^ [^; the fact that (in this case, at least for 
our smaller systems) the two peaks can even be displayed on the same scale is a reflection 
of the exceptionally delicate balance between the two entropy densities. 

Figure (^ allows one to see that fee is the thermodynamically preferred structure. This 
conclusion is expressed quantitatively in the results gathered in Table |1[ Our results at 
p = 0.7778 correct those of our earlier work |3^, as explained in section |III A| . They are 
in full accord with the results (both LS and IM-based) reported by Pronk and Frenkel 



3T| . The close correspondence between the results for N = 1728 and A^ = 5132 confirms 
that the former system is already representative of the thermodynamic limit. Table |I| 
also shows the results of our studies at the close-packed limit, using the hard-dodecahedron 



representation (Eq. |10a| ). Our results seem at variance with the IM-based result of Woodcock 



Q, even allowing for the large uncertainty attached to that result. They are close to those 



(based on LS) reported by Mau and Huse [|I^]. But the differences (for the smaller systems. 



particularly) appear to be statistically significant ||6^. Figure (^ gives an alternative view 
of these results. It utilizes the parameterization of the measured pressure difference between 
the two phases provided by Speedy [Q to determine the entropy difference, as a function of 



density, given the entropy difference at a chosen reference density; we have used the results 
of the present work at p = 0.7778. 



Table |T| shows the results of our studies in the constant pressure ensemble. The quantity 
of interest here is the difference between the gibbs free energy densities at the chosen pressure, 
which follows from the relevant distribution with the aid of Eq. (plsj) . In fact the gibbs free 
energy density difference Ag for a given pressure, and the entropy density difference As at a 
physical density that is the thermodynamic conjugate of that pressure for one of the phases, 
differ (in magnitude |6^) by terms that are seeond order in the pressure difference between 



the two phases. That pressure difference is extremely small [Q , as is the difference between 
the measured densities of the two structures (Table p|) . In these circumstances one would 
expect the magnitude of Ag to fall on the As plot in Figure (|]); and indeed, within the 
residual uncertainties, it does. 



V. DISCUSSION: REVIEW AND PROSPECTS 

In the work described here we have been concerned both with a system of long-standing 
interest -the hard sphere crystal- and a mei/ioc? -lattice-switch Monte Carlo- with potentially 
wide applicability. We divide our concluding discussion accordingly. 

The full agreement between the present work and that of [^ leaves little doubt that 



the equilibrium entropy difference between the two close-packed structures has finally been 
established securely and with high precision — at least at one density. Although a small 
discrepancy with respect to the results of |T^ remains, the accord of our close-packed limit 



results with those established using pressure difference measurements [Q suggests that the 
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curve in Fig. (j^) provides a relatively complete and trustworthy picture of the density- 
dependence. 

Notwithstanding the simplicity of the model, these results do have implications for 
experimentally-realizable systems. The immediate relevance to atomic systems is tenuous 
||68|| , but the model has been widely used to account for the behavior of assemblies of 'hard', 
'spherical' colloidal particles [^. Since the predicted entropy-density difference is so small 
there are potentially many ways (residual interactions between the spheres; polydispersity) 
in which the applicability of the theory may be compromised. But, of these, it seems that 
the most significant issues to be addressed are to do with scales -length and time. 

First, the lengths. In the experiments reported in ||6^ the colloidal particles have diam- 



eters of order 10^^ m and the samples comprise crystallites with linear dimensions of order 
10~^m. The number of particles in such a crystallite (A^ ~ 10^) is large compared to those 
in our simulation, which is (as we have seen) sufficient to allow us to deduce properties of the 
thermodynamic limit. But it is not large enough to guarantee that the behavior displayed 
will actually be that of the thermodynamic limit. To see this -and its principal implications- 
one needs to consider the stability of the perfect fee crystal with respect to hcp-tjpe stacking 
faults. Following reference [^ we may introduce a parameter a [^ measuring the proba- 
bility that a chosen close-packed plane sits within an fee environment as distinct from the 
hep environment. A simple argument (Appendix A) using the pseudo-spin parameterization 
of stacking patterns provided in |jl3] then yields the result 



« = -(! + tanh 
2 ' 



NiAs 



(24) 



where A^_l is the number of particles in a close packed layer and As (a function of p) is 
the fee-hep entropy difference per particle, as given in (and in the units of) Table |I|. The 
thermodynamic ideal [a = 1) is thus realized only to the extent that A^_l As is large compared 
to unity. For the length scales given above, A^_lAs ^ 1. The obvious implications are 



qualitatively consistent with the observations reported in |6^ which show a values (deduced 
from Bragg scattering intensities) ranging from 0.5 (signaling essentially random-hexagonal- 
close packing, rhep) through to a = 0.8. 

The observed spread in a values reflects -presumably- the issue of time scales. The 
smallness of the entropy difference (which supplies the kinetic driving force towards the 
equilibrium state) suggests that the equilibrium behavior will be observed only in samples 
which are grown sufficiently slowly and (or) given sufficient time for subsequent annealing 



71| . The results of ||6^ do indeed suggest a correlation between observed a value and 
the slowness of the growth process. Experiments done in microgravity [0, where growth 
processes are greatly accelerated, yield essentially randomly close-packed crystals. 

Now let us turn to the lattice-switch method. There are two questions here: One: does 
the method represent a significant advance with respect to existing methods? Two: is it 
generally applicable? 

The main alternative method (the benchmark against which others need to be assessed) 



is probably integration along a reference path, of which the work reported in |^ represents, 
to our knowledge, the most refined example. If one compares the two techniques (LS and 
IM) on the basis of precision-for-computational-buck there seems to be no clear winner 
in the hard-sphere studies to date: reference [^ reports calculations using both methods 



16 



that achieve comparable levels of precision on the basis of comparable computational time. 
But one should note that the entropy difference ultimately determined is some four orders of 
magnitude smaller |7^ than the separate entropies of the two phases, determined by IM. One 



can see this as a testimony to the care with which the recent IM studies have been carried out; 
or (as suggested in section |I[) as a strong indicator that another approach using an inter-phase 
path is called for. There are also two other counts — both somewhat subjective — on which 
we suggest that the LS approach is superior. First, it seems to us relatively illuminating (by 
comparison with IM) to read-off the result for a free energy difference directly from a figure 
like Fig. (^ which shows what it means. Secondly it also seems to us that LS wins in regard 
to the transparency of the uncertainties to be attached to its results. The LS error bounds 
represent purely statistical uncertainties associated with the measurement of the relative 
weights of two distribution-peaks. The IM error bounds have to aggregate the uncertainties 
associated with different stages of the integration process. 

As regards the second question, we expect that the method will, with appropriate ex- 
tensions, be widely applicable. The first extension must clearly be to accommodate soft 
potentials. The LS operation will then need gateway configurations in which the energies 
of the two structures (measured with respect to their ground-state energies — or indeed any 
fixed reference energy) are closely- matched [Q. The 'overlap order parameter' will need to 



be redefined accordingly. With no more than this degree of elaboration the method should 
be applicable immediately to investigate the widespread 'competition' between fee and hep 
ordering in the phase-behavior of the elements [^ . 



More generally, moving beyond the space of fcc-hcp structures, the choice of lattice-to- 
lattice mapping will require some thought. Mappings which preserve the relative positions 
of significant subsets of the particles (the analogues of the close-packed planes) are likely to 
be be optimal. The license to choose ones representation of the displacements (Section [11(JD 
may also prove useful. Simple transformations ||5^ will help if the mapping takes particles 
between environments in which the spectrum of single-particle displacements is significantly 
different. In such cases one might envisage using a MC-annealing procedure to refine the 
choice of representation. The use of normal coordinates has some advantages here -but pos- 
sibly not enough to offset the fact that the interaction potential is non-local when expressed 
in fourier space. 
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APPENDIX A: DISPLACEMENT ENTROPY VERSUS STACKING ENTROPY 

Consider a system of A^ hard spheres arranged in N\\ close-packed layers of N± particles. 
Following reference [Q one may conveniently index each of the close-packed layers with a 
pseudo-spin (Ising-like) variable a, where ctj = 1 signifies that layer i has an fee environment 
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(the two immediately adjacent layers are not aligned with one another) while ctj = — 1 implies 
an hep environment (adjacent planes are aligned with one another). The probability of a 
particular stacking sequence {a} (if these variables may properly be regarded as annealed) 
then satisfies 



In P{{a}\N, V) = S{N, V, {a}) + constant 



(Al) 



where S{N, V, {a}) measures the entropy associated with the configurations (displacements) 
consistent with the particular structure {a}. Following WM this entropy (we will refer to it 
here as 'displacement entropy') can usefully be written in the form of an expansion: 



S{N, V, {a}) = Nso + N^h^ai + N^J ^ a^a,- + . . . 



(A2) 



The expansion is effectively ordered in the range of the entropic inter-layer 'interac- 
tions': the ellipsis represents contributions from interactions (microscopically, displacement- 
displacement correlation functions) extending over more than 4 layers. The analysis of 



reference [|14| indicates that the series converges quickly, except close to melting. If we 
neglect the interaction terms altogether we may make the identification 



h 



and, from Eq. ( |Al|) , 
1 



[S{N, V, {a = +!})- SiN, V, {a = -1})] 



Asf. 



cc,hcp 



(A3) 



(^) 



A^ 



J2 ^({^}|^, V)ai = tanh [N±h] = tanh 



Wh 



NiAs 



fcc,hcp 



(A4) 



from which Eq. (|2^) follows. The correspondence with a ID paramagnet is clear. The famil- 
iar competition (between orientation energy and entropy) is played out here as a competition 
between displacement entropy and stacking entropy, with A^_l playing the role of an inverse 
temperature. 



TABLES 



mapping description effect m = Ai/N 

i (0, -t, +i) fee -^ hep 0.150(1) 

2 (0,2f,-2t) fee ^ hep 0.183(1) 

3 (0,3f,-3t) fee ^ fee 0.194(1) 

4 random-plane fee -^ hep 0.373(2) 

5 random-site fee -^ hep 0.820(3) 

TABLE I. The efficiency of different lattice mappings (for N = 1728 and p = 0.7778), as 

measured by the number of overlaps (per sphere) that they generate. Refer to the text for details. 
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p/p, 



cp 



N 



As (10-5 X k 



B) 



Method 



Ref. 



0.731 


512 


85 


0.736 


12000 


230 


0.736 


12096 


87 


0.739 


512 


90 


0.7778 


216 


132 


0.7778 


1728 


112 


0.7778 


1728 


113 


0.7778 


216 


133 


0.7778 


1728 


113 


0.7778 


5832 


110 


1.00 


12000 


260 


1.0 


512 


110 


1.0 


64 


91 


1.0 


216 


107 


1.0 


512 


119 


1.0 


1000 


113 


1.0 


216 


131 


1.0 


1728 


125 



(10) 
(100) 

(20) 
(4) 
(4) 
(4) 
(4) 
(3) 
(3) 
(3) 
(100) 

(20) 
(5) 
(4) 
(3) 
(4) 
(3) 
(3) 



SM 
IM 
IM 
LS 
LS 
LS 
IM 
LS 
LS 
LS 
IM 
SM 
LS 
LS 
LS 
LS 
LS 
LS 





m 
m 
ig 
11 
11 



PW 
PW 
PW 

m 

|14| 












PW 
PW 



TABLE II. The difference in the entropy densities of the fee and hep structures, 
As = As/^g f^ (Eq. |ll|); the associated uncertainties are in parenthesis. The results of the 



present work (PW) supercede those of reference |3C]. The results of reference |64| supercede those 
of reference ||in|. IM stands for integration method; SM is the lattice shear method of [14,|63[. 
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P* (D-3) 


Phcp 


Pfcc 


N 


A5 (10-5 X kBT) 


14.58 
14.58 


0.7776(1) 
0.7770(3) 


0.7775(1) 
0.7774(2) 


216 
1728 


-113 (4) 
-115 (6) 



TABLE III. The difference in the gibbs free energy densities of the fee and hep structures 
IS.g = ^gf(.(. fi(-„ (Eq. ^); the associated uncertainties are in parenthesis. P* gives the pressure 
I in units of ksT/D^. 
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FIGURES 






FIG. 1. Schematic representations of the different ways in which multicanonical samphng 
methods can be used to achieve inter-phase crossing. In the conventional approach (a) the samphng 
algorithm is biased so as to enhance the probability of the mixed phase states lying along a path 
(the heavy dark line) linking the two regions of configuration space. In the lattice-switch method 
(b) the bias is constructed so as to enhance the probability of the subsets of states (the white 
islands), within the single-phase regions, from which the switch operation (the large dashed arrow) 
will be accepted. 
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FIG. 2. Schematic representations of the two close-packed structures. The structures differ 
only in regard to the stacking pattern of the close-packed {x-y) planes which are of the form 
ABC ABC . . . for fee (upper) and ABAB . . . for hep (lower). The vector labeled t is instrumental 
in defining the LS operation, shown in Fig. (^. 
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FIG. 3. The LS transformation applied to the perfect- crystal configuration. The diagram 
shows 6 close-packed {x — y) layers. [The additional bracketed layer at the bottom is the periodic 
image of the layer at the top.] The circles show the boundaries of hard spheres located at the sites 
of the two close-packed structures. In this realization of the fee -^ hep lattice-switch, the top pair 
of planes are left unaltered, while the other pairs of planes are relocated by translations, specified 
by the vectors —t (white arrows) and t (black arrows). The vector t is identified in Fig. (^). 
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FIG. 4. The LS transformation applied to a 'typical' configuration. The crosses identify the 
'lattice sites'; the small circles locate the sphere centers in this configuration of displacements {u}- 
This configuration is realizable (gives no overlaps) in the fee structure; under the LS transfor- 
mation it is mapped onto an (unrealizable) hep configuration with four overlapping pairs of hard 
spheres (shown with dashed boundaries). Thus, for this configuration, the overlap order parameter 
M{{u]) = 4 (Eq. H). 
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FIG. 5. Illustration of the weight-generation process, for a system of A^ = 216 hard spheres. 
The points marked VS are the results of the first 3 iterations of the visited-states algorithm, 
initiated from an fee equilibrium state. The points marked TP emerge from one application of the 
transition probability method. The solid line shows a refined (usable) set of weights. 
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FIG. 6. Distribution of the separation d between adjacent close-packed planes in a system 
of 216 spheres at p = 0.7778, in the equilibrium hep and fee macrostates, and in the gateway 
[M = 0) macrostate. The separation is measured with respect to the equilibrium separation dg 



and is expressed in units of the sphere separation 5 f^\. 
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FIG. 7. (a) The mean value of the separation d between adjacent close-packed planes in a 
system of 216 spheres at p = 0.7778, for macrostates of different M. . The separation is measured 
with respect to the equilibrium separation d^ in units of 5 |^. Category-(i) planes (see Fig. ^) are 
translated together by LS; category-(ii) planes are translated through different amounts by LS. 



(b) The evolution with M of the c/a-ratio |55] in a constant-presswre ensemble (with the same 
parameters as (a)). The horizontal line marks the ideal-close-packed value. 
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FIG. 8. The probability distribution of the overlap order parameter per particle, m = M/N, 
for systems of three different A^ values at p = 0.7778. The lines provide gaussian guides to 
the eye; the statistical uncertanties on the data points are smaller than the symbol size. The 
entropy-difference is identified from the logarithm of the ratio of the integrated weights of the two 
peaks. The hep peak for the largest system is not visible on this scale. 



29 



0.0014 



0.0012 - 



0.001 



CO 

< 



0.0008 



0.0006 



[] 



DMau & Huse(1999) 
OPronk&Frenkel(1999) 
• PW 



0.7 



0.8 



P/P 



0.9 



CP 



FIG. 9. The difference in the entropy densities of the fee and hep structures, As = As^^^ ^g„ 
(Eq. [l^), as a function of reduced density p. The data points are as given in Table ||. The solid 



line is the result of an integration of the pressures of the phases [66|. Note that this line passes 
through our result at p = 0.7778 by construction. 
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